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A unified method to permute vector-stored Upper-triangular-Diagonal factorized 
covariance (UD) and vector-stored upper-triangular Square-Root Information (SRI) arrays 
is presented . The method involves cyclic permutation of the rows and columns of the 
arrays and retriangularization with fast (slow) Givens rotations (reflections). Minimal 
computation is performed , and a one-dimensional scratch array is required. To make the 
method efficient for large arrays on a virtual memory machine , computations are arranged 
so as to avoid expensive paging faults. This method is potentially important for processing 
large volumes of radio metric data in the DSN. 


I. Introduction 

In the reduction of observational data involving simultane- 
ous least-squares estimation of many parameters, it is desirable 
to examine different modeling scenarios without reprocessing 
the data through the filter. Often the analyst will estimate 
a full contingent of system parameters in the initial filter pass. 
A reduced state estimate, reflecting a different modeling sce- 
nario, can then be obtained from this initial pass. Furthermore, 
the sensitivity of the reduced state to the excluded parameters 
is readily obtainable. Application of a priori information for 
the excluded parameters to the sensitivity and subsequent aug- 
mentation of the reduced state by these perturbations then 
yields a more conservative assessment of the reduced state 
errors. This type of analysis is particularly useful when the 


reduced state is sensitive to the excluded parameters that are 
insensitive to data [1] . 

Upper-triangular data structures, such as those encountered 
in an Upper-triangular-Diagonal factorized covariance filter 
(UD filter) or in an upper-triangular Square-Root Information 
filter (SRI filter), are conducive to computing these reduced 
state estimates. In a UD filter, it is necessary to permute the 
rows of the UD array in order to compute the reduced esti- 
mate. The rows corresponding to the parameters that are to be 
excluded must be moved below those corresponding to the 
parameters that are to be included in the reduced estimate. In 
an SRI filter, it is necessary to permute the columns of the 
SRI array in order to compute the reduced estimate. The 
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columns corresponding to the parameters that are to be 
excluded must be moved to the right of the columns corre- 
sponding to those that are to be included in the reduced 
estimate. 1 

A brute force method to accomplish these permutations 
reorders the rows (columns) of the UD (SRI) array into a 
two-dimensional work array. To retriangularize this two- 
dimensional array, one can then post(pre)-multiply it by an 
implicitly defined orthogonal transformation composed of 
Householder reflections. 2 This method is sound and numeri- 
cally stable due to the orthogonality of the transformation. 
Furthermore, the commercially available Estimation Subrou- 
tine Library (ESL) has standard routines to conveniently per- 
form these permutations and subsequent retriangularizations 
[2] . However, this method requires a two-dimensional work 
array to store the permuted UD (SRI) array. For large-scale 
systems, such as those encountered for Global Positioning 
System (GPS) studies, it is not at all uncommon to estimate on 
the order of 1000 parameters. This translates into four million 
extra bytes of storage. Additionally, over 3000 parameters are 
routinely solved when DSN VLBI data are processed with 
MASTERFIT software [3] . Moreover, for systems of this 
order, retriangularization is computationally expensive. 

In 1986, Bierman (personal communication) suggested that 
the necessary reordering of UD arrays be performed as a series 
of pairwise permutations in such a way that the upper- 
triangular structure of the UD array is always maintained. Each 
pairwise permutation could be performed as a suboptimal 
noise-free measurement update: 

Update = (/ - GH) P(I - GH) 1 + RGG* 

where P update is then the permuted covariance of the covari- 
ance P. The noise R is set equal to zero, and the gain G and 
design H are chosen so that (/ - GH) is a permutation opera- 
tor. Such a permutation operator that permutes parameters i 
and / can be expressed as 

!-GH = /-(e.-e.)(e.-e./ 


J S. C. Wu et al., Oasis Mathematical Description , JPL Publication 
D-3139 (internal document), Jet Propulsion Laboratory, Pasadena, 
California, April 1, 1986. 

2 Before this retriangularization can be performed on the UD array, it 

is of course necessary to scale the columns by the square root of the 

corresponding diagonal. These square roots can be avoided by per- 

forming a modified weighted Gram-Schmidt UD triangularization 

algorithm. 


where e { is a column vector whose elements are zero except 
element /, which is 1. Trivially H - (e- - e^Y and G = H x . In 
practice this computation is implemented with an optimal 
measurement update using the Bierman UD measurement 
update algorithm and a rank-1 update to include the effect 
of the suboptimality of the gain. This method eliminated the 
need for a two-dimensional scratch array. Software was even- 
tually written that exploited the structure of the suboptimal 
gain G and design H. This method still required several scratch 
arrays and a multitude of pairwise permutations. 3 

In 1987, Wolff (personal communication) suggested that 
the same pairwise permutations could be performed on an SRI 
array with retriangularization accomplished after each column 
exchange with slow Givens reflections. 4 This method also 
eliminated the need for a two-dimensional work array. How- 
ever, in the retriangularization process, elements not consistent 
with the upper-triangular data structure of the SRI array 
were created. Each inconsistent element required another slow 
Givens reflection. Subsequently, this method required on the 
order of n 2 slow Givens reflections, where n is the column dis- 
tance between the two parameters of the exchange. Further- 
more, this method required two explicit scratch arrays. One of 
the scratch arrays stored the right column of the exchange 
while the other was used to store elements that were created in 
the retriangularization process. 5 

Also in 1987, Pombra suggested a technique that eliminated 
the need to perform a multitude of these pairwise permuta- 
tions in reordering a UD array. 6 His technique was to insert 
into the covariance an artificial parameter with zero variance 
and zero correlation with the other parameters. This is easily 
accomplished by inserting into the UD array a row and column 
of zeros. A parameter could then be moved directly into its 
proper place by performing a pairwise permutation of this 
artificial parameter and the parameter that is desired to reside 
in its place. After the permutation is performed, this artificial 
parameter can then be deleted from the UD array. This is also 
easily accomplished by removing from the UD array the appro- 
priate row and column of zeros. This technique drastically 


3 P. J. Wolff, “UD Permutations via ESL Subroutine UCON,” JPL IOM 
314.5-1040 (internal document), Jet Propulsion Laboratory, Pasa- 
dena, California, August 29, 1986. 

4 Since the SRI array is not unique, either Givens reflections or rota- 
tions can be used. 

5 P. J. Wolff, “Permuting R Matrices in Place,” IOM 314.5-1091 (inter- 
nal document), Jet Propulsion Laboratory, Pasadena, California, 
March 10, 1987. 

6 S. A. Pombra, “Computationally Fast Version of UDEDIT,” IOM 
335.1-87-160 (internal document), Jet Propulsion Laboratory, Pasa- 
dena, California, July 7, 1987. 
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reduced the number of permutations necessary to reorder the 
array. However, suboptimal noise-free measurement updates 
were still used to perform the pairwise permutations. 

It was recognized that the insertion of a row and column of 
zeros into the UD array, together with the subsequent row 
exchange between the artificial parameter and the parameter 
of interest, created an array similar to a Morf-Kailath prearray 
with zero noise variance [4] . Retriangularization of this array 
into an upper-triangular Morf-Kailath postarray is easily 
accomplished with fast Givens rotations. 7 This method re- 
quires a one-dimensional scratch vector which manifests itself 
as an augmentation of the UD array. An extra column of the 
UD array permits the insertion of a row and column of zeros. 
Furthermore, the method requires on the order of 7T row fast 
Givens rotations, where fr row is the row distance between the 
parameter and the desired location of the parameter. Due to 
the nature of the fast Givens, square root computations are 
never required. 


II. Row Permutations of UD Arrays 

To demonstrate the technique on UD arrays, let the factori- 
zation of a covariance P with parameters Ul, U2, U3, U4, U5 
be 


where 


P = UDU f 


wl l 
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u 2 2 

U ~3 

" 2 4 

“ 2 S 
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0 
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0 
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( 1 ) 


( 2 ) 


This technique can also be adopted to reordering SRI 
arrays. In a manner similar to the Wolff method, retriangulari- 
zation is accomplished with slow Givens reflections. Unlike the 
Wolff method, however, no new elements are created. Hence, 
this method requires on the order of only n co \ umn slow Givens 
reflections, where Column is the column distance between the 
parameter and the desired location of the parameter. 

Since the artificial parameter is only a convenient place 
holder, an alternative method could be constructed that does 
not incorporate the insertion/deletion step. Similar results can 
be achieved by cyclically permuting the rows and columns 
of the arrays. Retriangularization is then achieved as before. 
Furthermore, since the size of the array does not change, it is 
not necessary to augment the array with an extra column. 
However, a one -dimensional scratch vector is required. 


and 


D = 


0 

0 


0 

d 2 

0 


0 0 
0 0 


0 0 0 

0 0 0 

d 3 0 0 

0 d 4 0 
0 0 d 5 


If we designate a weighted inner product as 


(3) 


Cyclic column permutations and retriangularization of SRI 
arrays have already appeared in the literature [5] . In fact, the 
discussed method to cyclically permute columns of an SRI 
array is similar to the left circular shift method used in the 
UNPACK routine SCHEX [6]. 8 Until now, however, these 
results have not been extended to cyclic row permutations and 
retriangularization of UD arrays. The permutations are per- 
formed in such a way that the only difference between per- 
muting UD and SRI arrays lies in the retriangularization 
process. 


7 Since the UD array is unique, fast Givens reflections cannot be used. 

8 SCHEX requires two one-dimensional scratch arrays to store the 
transformation pairs and does not store the upper-triangular array as a 
vector. 


* = Tj x i d i y i (4) 


then the covariance can be written in terms of its factorized 
elements: 


m 1 • Ml 
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m5 • m4 

m5 • m5 


(5) 
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Rows 2 and 5 of Eq. (7) can then be exchanged: 


As usual, let ui i = 1 to obtain a unique factorization. With this 
convention, the diagonals of the D array can now be stored as 
the diagonals of the U array. This new array is referred to as 
the UD array: 9 
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Say it is desired to reorder the parameters as U 1, U 3, U 4, 
U 2, US. The first step is to insert an artificial parameter, with 
zero variance and zero correlation with the other parameters, 
before parameter US. This is accomplished by inserting into 
the UD array a row of zeros above row 5 and a column of zeros 
left of column 5: 


In practice, however, elements d 2 , m2 3 , and m2 4 would not be 
moved, since a scratch array would be needed to store these 
elements. As will soon be evident, the second row can itself be 
used as a scratch array. The array in Eq. (9) implicitly corre- 
sponds to a nontriangular U array 
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The corresponding covariance is 
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9 The upper-triangular elements of either a UD or an SRI array are 
stored contiguously in computer memory as a vector; the zero ele- 
ments below the diagonal are obviously not needed to represent the 
array. 
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and a covariance P 

p= mu 1 
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Except for the artificial parameter, this is the desired permuted 
covariance. The covariance factorization in Eq. (12) can be 
thought of as the product of an upper-triangular array and its 
transpose: 


P = V s/ffiVsJB)* 


(13) 


where 
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It is desirable to construct an orthogonal transformation T 
such that U >JBt is an upper-triangular array with zeros in the 
second row and column: 

Us/St = 


Us/Bt = 


ul 1 y/d 1 

0 

0 

0 

0 

0 


0 x 
0 0 
0 * 
0 0 
0 0 
0 0 


0 

X 

X 

0 


X 

0 

X 

X 

X 


0 0 


Ml sV^ 5 

0 

«3 5 v^ 5 

u4 s y/d s 

u2 s s/d s 

u5.\/d 


( 16 ) 


A portion of the array in Eq. (15) is similar to a Morf-Kailath 
prearray with zero noise variance. Reduction of this prearray 
to the postarray in Eq. (16) is easily accomplished with three 
fast Givens rotations. As a result of the zero noise variance, the 
first fast Givens T 0 is simply a permutation of columns 2 
and 5: 
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(17) 
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This is the required step to produce the zeros in the second 
row and column. The next fast Givens T x operates on columns 
3 and 5 in such a way that the array element containing 
u 2 3 \47 3 becomes zero. The final fast Givens T 2 operates on 
columns 4 and 5 such that the array element containing 
u2^\/d 4 becomes zero. Each of these Givens transformations 
maps the elements in each row of the columns on which they 
operate into a linear combination of the same. Hence, ele- 
ments not consistent with an upper-triangular structure are 
never created; a linear combination of zeros is still zero. For 
example, 7\ operates on Eq. (17) as follows: 


Likewise, T 2 operates on Eq. (18) as follows: 
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Trivially, since the transformations 7} are orthogonal, the 
covariance Pin Eqs. (12) and (13) remains unchanged: 


F= Us/Btjj^Us/EtjjJ 
= U BtP (20) 


The importance of the zeros in the second row and column of 
Eq. (16) is now readily apparent. Since the insertion of a row 
and column of zeros into the UD array in Eq. (7) did not 
affect the individual elements of the covariance in Eq. (8), 
deletion of the second row and column of zeros from Eq, (16) 
will have a similar null effect on the individual elements of the 
covariance P. 
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III. Column Permutations of SRI Arrays 

To demonstrate the technique on SRI arrays, represent an 
augmented SRI array R with parameter order R 1, R2 , R3, 
R4, R 5 as 
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Columns 2 and 5 of Eq. (23) can then be exchanged. This is 
just a reordering of the design equation 
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( 21 ) 

where Z is the residual vector. The following identity suggests 
the justification of incorporating an artificial parameter into 
an SRI array: 
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If there are no intentions of computing the inverse of the SRI 
array, it is permissible to allow e = 0. Hence, the insertion and 
subsequent deletion of an artificial row and column of zeros 
into an SRI array does not affect the information content of 
the array. 


Say it is desired to change the parameter order to/?l,/?3, 
R 4, R 2, R 5. The first step is to insert a row and column of 
zeros into the (R \Z) array before the parameter R 5: 


It is desirab le to co nstruct an orthogonal transformation T' 0 
such that T 0 (R\Z ) has only zeros in the second row and col- 
umn. Such a transformation is just a permutation of rows 2 
and 5: 
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In practice, however, elements r 3 2 and r4 2 would not be 
moved, since a scratch array would be needed to store these 
elements. As is evident, the second row can itself be used as a 
scratch array. It is now a simple matter of retriangularizing 
Eq. (25) with two slow Givens. The first Givens T x operates 
on rows 3 and 5 in such a way that the array element contain- 
ing r 3 2 becomes zero. The second Givens T 2 operates on 
rows 4 and 5 in such a way that the array element containing 
r 4 2 becomes zero. Each of these Givens transformations maps 
the elements in each column of the rows on which they oper- 
ate into a linear combination of the same. Elements not con- 
sistent with an upper-triangular structure are never created. 
For example, T l operates on Eq. (25) as follows: 

1 0 0 0 0 (f 

0 1 0 0 0 0 

0 0 c 0 s 0 

0 0 0 1 0 0 

0 0 ±s 0 +c 0 

0 0 0 0 0 1 


Likewise, T 2 operates on Eq. (26) as follows: 
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(27) 



’ rl l 

0 

' 3 1 

r4 l 

r2 l 

r5 i 

Z i 


0 

0 

0 

0 

0 

0 

0 


0 

0 

a: 

X 

X 

X 

X 


0 

0 

0 

r4 4 

0 

r5 4 

Z 4 


0 

0 

0 

X 

X 

X 

X 


0 

0 

0 

0 

0 

r5 s 



( 26 ) 


From Eq. (22) it is a trivial matter that the second row and 
column of zeros of Eq. (27) can then be deleted. 


IV. Cyclic Permutations 

Since the artificial parameter is only a convenient place 
holder, an alternative method can be constructed that does not 
incorporate the insertion/deletion step. By cyclically per- 
muting the rows and columns of either the UD or SRI arrays, 
a data structure similar to Eq. (17) or Eq. (25), respectively, 
can be achieved. For example, cyclically left permuting col- 
umns 2 through 4 of Eq. (21) yields the following: 
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Cyclically upward permuting 
then yields 
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rows 2 through 4 of Eq. 

(28) 
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V. Implementation 

When dealing with large arrays on a virtual memory ma- 
chine, computations should be arranged to minimize expensive 
page faulting. The Givens sweeps should access contiguous 
storage locations. The usual column ordering of a vector-stored 
upper-triangular data structure with n rows and n columns is 
as follows: 


column 


®(1) 

s(2) 

5(4) 

0 

»(3) 

5(5) 

0 

0 

5(6) 

0 

0 

0 





This is the desired data structure of Eq. (25) without the arti- 
ficial parameter. Slow Givens sweeps across rows 2,4 and 
rows 3,4 of Eq. (29) can then be applied to zero out the array 
elements containing r 3 2 and r 4 2 , respectively. Since the size 
of the array never changes, this method would not require that 
the array be augmented with an extra column. However, one 
explicit scratch array is required to hold that part of Eq. (29) 
not consistent with an upper-triangular data structure. 


(31) 

Element S column (/, /) is easily referenced as s [(/ 2 - j/2) +/]. 
For UD arrays, the Givens sweeps down the columns using the 
data storage in Eq. (31) would require minimal paging. For 
SRI arrays, the Givens sweeps across the rows with the data 
storage in Eq. (31) would be paging intensive. Hence, for 
SRI arrays, it proves to be more efficient to store the upper- 
triangular array with unnatural row ordering as follows: 


A similar approach can also be taken toward permuting 
UD arrays. Cyclic permutations of the rows and columns 
would yield the data structure in Eq. (17) without the artifi- 
cial parameter: 
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(30) 

Fast Givens sweeps down columns 2,4 and columns 3,4 of 
Eq. (30) can then be applied to zero out the array elements 
containing w 2 3 V ^3 an ^ respectively. 
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Element S roy/ (z, / ) can be referenced as s[ j + ( / 2 - i/2) 
+ (/-!)(«- i)]. 10 Since filtering software generally stores 


10 ^rowO’> /) “ ^columnO’ / )J this is just a relabeling of the elements. 
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upper-triangular arrays according to Eq. (31), it would be 
necessary to relabel from S column to S row before the permu- 
tations and then back to S column afterward. 

It is equally easy to perform cyclic permutations with 
either Eq. (31) or Eq. (32). In the former, movement is along 
the columns; in the latter, movement is across the rows. The 
only paging-intensive step in cyclically permuting the rows and 
columns with data storage S column (^ row ) * s * n stor i n 8 the 
upper row (left column) into a scratch array. 

VI. Performance 

To assess the permutation methods on UD arrays, an upper- 
triangular array was constructed with 499 parameters and an 
estimate column. The following permutations were then per- 
formed on a VAX 11/780. First, the parameters were per- 
muted randomly. In practice, of course, this would never hap- 
pen. It is more realistic to move blocks of parameters to the 
bottom of the array. Therefore, the following systematic per- 
mutations were performed. A block of 10 parameters starting 
in rows 1, 167, and 334 were moved to the bottom of the 
array. These permutations were then repeated with a block of 
100 parameters. CPU time and page faults were accumulated 
for three different methods. First, the ESL routine HHPOST 
was used to apply a post-Householder to a row-reordered 
column-scaled UD array. Second, pairwise permutations were 
performed with suboptimal noise-free measurement updates. 
This is equivalent to the ESL routine U2U. Third, the permu- 
tations were performed with the discussed method of cyclic 
permutations with fast Givens sweeps. The results are pre- 
sented in Table 1. 

The ESL method with HHPOST performs the same regard- 
less of the number and type of permutations performed. The 
performance of the other methods depends on the number and 
type of permutations. The cyclic permutation method with 
fast Givens sweeps is shown to be superior to the ESL method 
in U2U. Furthermore, in the cyclic permutation method, as 
the starting row moves down there are fewer elements to 
retriangularize. 


The ESL has standard routines to perform permutations on 
SRI arrays. The routine R2A copies and reorders the columns 
of an upper-triangular array into a two-dimensional work array. 
A Householder transformation can then be applied with rou- 
tine TDHHT to retriangularize back to an upper-triangular 
array. 

To assess the permutation methods on SRI arrays, ESL rou- 
tines R2A/TDHHT, the discussed method of cyclic permuta- 
tions with slow Givens sweeps using data storage S co lumn , and 
the same with data storage S row were used to permute an 
upper-triangular array with 499 parameters and a residual 
column. CPU time and page faults were accumulated for per- 
muting, retriangularizing, and any necessary relabeling. The 
results are presented in Table 2. 

The ESL method performs the same regardless of the num- 
ber and type of permutations performed. The performance of 
the cyclic permutations with slow Givens sweeps depends on 
the number and type of permutations. As the starting column 
is moved to the right, not only are there fewer elements to 
retriangularize, but the lengths of Givens sweeps also become 
smaller, and fewer computations are required. Since the VAX 
11/780 is a virtual memory machine, the method with data 
storage S row is generally superior to the method with data 
storage S column . For fixed memory machines or small SRI 
arrays, the method with data storage S column should be 
superior, since no relabeling or unlabeling is required. 


VII. Conclusion 

The most efficient method to permute UD (SRI) arrays has 
been shown to be cyclic permutation with fast (slow) Givens 
sweeps for retriangularization. Since the only difference is in 
the retriangularization process, it has been possible to combine 
the methods into one subroutine. This method has been incor- 
porated into the OASIS software [7] and should prove to be 
advantageous for GPS, TOPEX, and DSN high-earth-orbiter 
data processing and analysis. 
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Table 1. CPU and paging faults in permuting UD arrays 



CPU time; thousands of page faults 

Permutations 

HHPOST 

ESL U2U 

Cyclic permutations 
and fast Givens 


0:33:06 

1:06:31 

0:17:04 

Random 

1025 

2051 

612 

Rows 

0:32:43 

0:48:55 

0:01:24 

1-10 

948 

1823 

40 

Rows 

0:31:49 

0:40:00 

0:01:14 

167-176 

970 

1338 

35 

Rows 

0:31:12 

0:24:26 

0:00:45 

334-343 

970 

728 

21 

Rows 

0:32:44 

0:59:25 

0:11:23 

1-100 

968 

1985 

363 

Rows 

0:31:32 

0:48:47 

0:08:46 

167-266 

946 

1433 

275 

Rows 

0:30:59 

0:27:29 

0:02:20 

334-433 

973 

747 

52 

Table 2. 

CPU and paging faults in permuting SRI arrays 


CPU time; thousands of page faults 

Permutations 

Cyclic permutations 

Cyclic permutations 

ESL and slow Givens 

and slow Givens 


TDHHT using data storage 
^column 

using data storage 
c 

row 


0:22:28 

0:41:09 

0:20:07 

Random 

729 

3253 

575 

Columns 

0:20:29 

0:01:36 

0:01:38 

1-10 

717 

48 

55 

Columns 

0:21:13 

0:00:55 

0:01:03 

167-176 

723 

31 

41 

Columns 

0:21:53 

0:00:20 

0:00:39 

334-343 

719 

11 

27 

Columns 

0:21:25 

0:13:11 

0:11:10 

1-100 

727 

424 

411 

Columns 

0:21:15 

0:06:14 

0:05:14 

167-266 

714 

229 

205 

Columns 

0:21:03 

0:00:51 

0:01:43 

334-433 

723 

52 

104 


